% This script produces Figures 1 and 2 from the paper.
% 07/2024: Colours have been adjusted to stand out under grey scale. 

clear
load('recalibrated_data_with_eval.mat')

% 1) Expected vs. actual wealth (Figure 1)
% This plots holdings of total and illiquid wealth over the life cycle,
% contrasting the actual holdings (_naif) with the holdings expected by a
% naive agent in the beginning of the life cycle (_TC).
% The variable capturing total wealth adds liquid and illiquid wealth
% holdings, net of current income flow. It is assumed that the agent holds
% 0.50 of their monthly income as part of their liquid wealth. 

age_ = mydata.age;
Zpath_naif = mydata.avgZ500_naif*1.905;
totalpath_naif = (mydata.avgX500_naif + mydata.avgZ500_naif - mydata.avgY + mydata.avgY/24)*1.905;
Zpath_TC = mydata.avgZ500_TC*1.905;
totalpath_TC = (mydata.avgX500_TC + mydata.avgZ500_TC - mydata.avgY + mydata.avgY/24)*1.905;
zpen_eval = mydata.zpen_eval;


figure
plot(age_(1:end-10),totalpath_naif(1:end-10),'k-', 'LineWidth',1.0)
hold on
plot(age_(1:end-10),Zpath_naif(1:end-10),'kx', 'LineWidth',1.0)
hold on
plot(age_(1:end-10),totalpath_TC(1:end-10),'r--','LineWidth',1.25)
hold on
plot(age_(1:end-10),Zpath_TC(1:end-10),'rx','LineWidth',1.25)
grid on
title(['Wealth accumulation over the life cycle ($ \beta = 0.503, \delta = 0.988, \theta = 1.105, R^Z = 5\%$) '],'interpreter','latex')
xlabel('Age')
ylabel('Wealth (2018 USD)')
legend({['Total wealth - actual'] ,['Illiquid wealth - actual'],['Total wealth - expected'] ,['Illiquid wealth - expected']},'interpreter','latex','location','northwest')


% 2) Savings paths, for different Rz (Figure 2)

% The following calculates a measure of the agent's forecasting error, as
% a function of the interest rate on the illiquid asset Z. It compares actual mean
% (positive) retirement savings to the naive agent's prediction. 
THETAmistake = [];

for THETA = mydata.THETA_grid
    suffix = strcat('avgZ',num2str(THETA*10000),'_naif');
    avgZ_naif = mydata.(suffix);
    suffix = strcat('avgZ',num2str(THETA*10000),'_TC');
    avgZ_TC = mydata.(suffix);
    
    Z_save_naif = zeros(1,70);
    Z_save_TC = zeros(1,70);
    
    % Back out period-by-period savings, taking into account withdrawal penalties.
    for t = 1:70
        if avgZ_naif(t+1)/(1+THETA) >= avgZ_naif(t)
            Z_save_naif(t) = avgZ_naif(t+1)/(1+THETA) - avgZ_naif(t);
        else
            Z_save_naif(t) = (avgZ_naif(t+1)/(1+THETA) - avgZ_naif(t))*(1-zpen_eval(t));
        end
        
        if avgZ_TC(t+1)/(1+THETA) >= avgZ_TC(t)
            Z_save_TC(t) = avgZ_TC(t+1)/(1+THETA) - avgZ_TC(t);
        else
            Z_save_TC(t) = (avgZ_TC(t+1)/(1+THETA) - avgZ_TC(t))*(1-zpen_eval(t));
        end
    end
    
    suffix = strcat('Zsave',num2str(THETA*10000),'_naif');
    workspace.(suffix) = Z_save_naif;
    suffix = strcat('Zsave',num2str(THETA*10000),'_TC');
    workspace.(suffix) = Z_save_TC;
    THETAmistake = [THETAmistake mean(Z_save_TC(Z_save_TC>0)) - mean(Z_save_naif(Z_save_naif>0))];
end


% The following produces Figure 2. 
figure
plot(age_(1:end-20),workspace.Zsave500_naif(1:end-19),'k-', 'LineWidth',1.25)
hold on
plot(age_(1:end-20),workspace.Zsave500_TC(1:end-19),'r-', 'LineWidth',1.25)
hold on
plot(age_(1:end-20),workspace.Zsave400_naif(1:end-19),'k--', 'LineWidth',1.25)
hold on
plot(age_(1:end-20),workspace.Zsave400_TC(1:end-19),'r--', 'LineWidth',1.25)
grid on
title(['Saving over the life cycle ($ \beta = 0.503, \delta = 0.988, \theta = 1.105$) '],'interpreter','latex')
xlabel('Age')
ylabel('Savings (2018 USD)')
ylim([-3e4 2.5e4])
legend({['Retirement saving - actual ($R^Z = 5\%$)'] ,['Retirement saving - expected ($R^Z = 5\%$)'],['Retirement saving - actual ($R^Z = 4\%$)'] ,['Retirement saving - expected ($R^Z = 4\%$)']},'interpreter','latex','location','southwest')

